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ABSTRACT 



A thorough search for large scale anisotropies in the distribution of arrival directions 
of cosmic rays detected above 10 18 eV at the Pierre Auger Observatory is presented. 
This search is performed as a function of both declination and right ascension in sev- 
eral energy ranges above 10 18 eV, and reported in terms of dipolar and quadrupolar 
coefficients. Within the systematic uncertainties, no significant deviation from isotropy 
is revealed. Assuming that any cosmic ray anisotropy is dominated by dipole and 
quadrupole moments in this energy range, upper limits on their amplitudes are de- 
rived. These upper limits allow us to challenge an origin of cosmic rays above 10 18 eV 
from stationary galactic sources densely distributed in the galactic disk and emitting 
predominantly light particles in all directions. 

Subject headings: astroparticle physics; cosmic rays 



1. Introduction 



Establishing at which energy the intensity of extragalactic cosmic rays starts to dominate the 
intensity of galactic ones would constitute an important step forward to provide further understand- 
ing on the origin of Ultra-High Energy Cosmic Rays (UHECRs) . A time honored picture is that the 



ankle, a hardening of the energy spec t rum located at ~ 4 EeV (ILinslev 19631 ; lLawrence et al. 1991 



Nagano et al. 19921 : iBird et al. 19931 : lAueer Collaboration 201031 ) (where 1 EeV = 10 18 eV), is the 
featu re in the ener gy spectrum marking the transition between galactic and extragalactic UHE- 
CRs (jLinslev 19631 ) . As a natural signature of the escape of cosmic rays from the Galaxy, large 
scale anisotropies in the distribution of arrival directions could be detected at energies below this 
spectral feature. Both the amplitude and the shape of such patterns are uncertain, as they depend 
on the model adopted to describe the regular and turbulent components of the galactic magnetic 
field, the charges of the cosmic rays, and the assumed distribution of sources in space and time. 
For cosmic rays mostly heavy and originating from st ationary sources loca t ed in the galactic d isk, 
some estimates based on diffusion and d rift motions (Ptuskin et al. 19931: Candia et al. 2003 ) as 
well as direct integration of trajectories ( Zirakashvili et al. 19981 : Giacinti et al. 20121 ) show that 
dipolar anisotropies at the level of a few percent could be imprinted in the energy range just below 
the ankle energy. Even larger amplitudes could result in the case of light prim aries, unless sources 
are strongly intermit tent and pure diffusion motions hold up to EeV energies (jCalvez et al. 20101 : 
Eichler & Pohl 2011 
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If UHECRs above 1 EeV have a predomina nt extragalactic origin (jHillas 19671 ; iBlumenthal 197Q ; 
Berezinskv et al. 2000 ; iBerezinskv et al. 2004 ). their angular distribution is expected to be isotropic 
to a high level. But, even for isotropic extragalactic cosmic rays, the translational motion of 
the Galaxy relative to a possibly stationary extragal actic cosmic ray rest fram e can produce a 
dipole in a similar way to the Compton- Getting effect (jCompton k, Getting 19351) which has been 
measured with cosmic r ays of much lowe r energy at the solar time scale (ICutler &, Groom 19861 : 
Amenomori et al. 20061 : lAbdo et al. 20091 : lAglietta et al. 20091 : lAbbasi et al. 2010a! ) as a result of 
the Earth motion relative to the frame in which the cosmic rays have no bulk motion. More- 
over, the rotation of the Galaxy can also produce anisotropy by virtue of moving magnetic fields, 
as cosmic rays travelling through far away regions of the Galaxy experience an electric force du e 
to the relative motion of the system in which the field is purely magnetic (jHarari et al. 20101 ). 
The large scale structure of the galactic magnetic field is expected to transform even a simple 
Compton-G etting dipole into a more complex anisotropy at Earth, described by higher order 
multipoles (jHarari et al. 20101 ). A quantitative estimate of the imprinted pattern would require 
knowledge of the global structure of the galactic magnetic field and the charges of the particles, 
as well as the frame in which extragalactic cosmic rays have no bulk motion. If, for instance, 
the frame in which the UHECR distribution is isotropic coincides with the cosmic microwave 
backg round rest frame, the amp litude of the simple Compton-Getting dipole would be about 
0.6% (jKachelriess &: Serpico 20061 ) . The same order of magnitude is expected if UHECRs have 
no bulk motion with respect to the local group of galaxies. 

The large scale distribution of arrival directions of UHECRs as a function of the energy is thus 
one important observable to provide key elements for understanding their origin in the EeV energy 
range. Using the large amount of data collected by the Surface Detector (SD) array of the Pierre 
Auger Observatory, results of first harmonic analyses of the right asc ension distribution performe d 
in different energy ranges above 0.25 EeV were recently reported (lAuger Collaboration 201 lal ). 
Upper limits on the dipole component in the equatorial plane were derived, being below 2% at 
99% C.L. for EeV energies and providing the most stringent bounds ever obtained. These analyses 
benefit from the almost uniform directional exposure in right ascension of the SD array of the 
Pierre Auger Observatory which is due to the Earth rotation, and they constitute a powerful tool 
for picking up any dipolar modulation in this coordinate. However, since this technique is not 
sensitive to a dipolar component along the Earth rotation axis, we aim in the present report at 
estimating not only the dipole component in the right ascension distribution but also the component 
along the Earth rotation axis. More generally, we present a comprehensive search in all directions 
for any dipole or quadrupole patterns significantly standing out above the background noise. 

Searching for anisotropies with relative amplitudes down to the percent level requires the 
control of the exposure of the experiment at even greater accuracy. Spurious modulations in the 
right ascension distribution are induced by the variations of the effective size of the SD array 
with time and by t he variations of the counting rate of events due to the changes of atmospheric 
conditions. In Ref. (lAuger Collaboration 201 lal ). we showed in a quantitative way that such effects 
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can be properly accounted for by making use of the instantaneous status of the SD array provided 
each second by the monitoring system, and by converting the observed signals in actual atmospheric 
conditions into the ones that would have been measured at some given reference atmospheric 
conditions. Searching for anisotropies explicitly in declination requires the control of additional 
systematic errors affecting both the directional exposure of the Observatory and the counting rate 
of events in local angles. Each of these additional effects are carefully presented in sections [3] and [H 

After correcting for the experimental effects, searches for large scale patterns above 1 EeV 
are presented in section [5j Additional cross-checks against eventual systematic errors affecting 
the results obtained in section [5] are presented in section [6j Resulting upper limits on dipole and 
quadrupole amplitudes are presented and discussed in section while a final summary is given in 
section El Some further technical aspects are detailed in the appendices. 



2. The Pierre Auger Observatory and the data set 



The Pierre Auger Observatory (| Auger Collaboration 20041 ). located in Malargiie, Argentina, at 
mean latitude 35.2° S, mean longitude 69.5° W and mean altitude 1400 meters above sea level, has 
been designed to collect UHECRs with unprecedented statistics. It exploits two available techniques 
to detect extensive air showers initiated by cosmic ray interactions in the atmosphere : a surface 
detector array and a fluorescence detector. The SD array consists of 1660 water-Cherenkov detectors 
laid out over about 3000 km 2 on a triangular grid with 1.5 km spacing. These water-Cherenkov 
detectors are sensitive to the light emitted in their volume by the secondary particles of the showers, 
and provide a lateral sampling of the showers reaching the ground level. At the perimeter of this 
array, the atmosphere is overlooked on dark nights by 27 optical telescopes grouped in 5 buildings. 
These telescopes record the number of secondary charged particles in the air shower as a function 
of depth in the atmosphere by measuring the amount of nitrogen fluorescence caused by those 
particles along the track of the shower. 

The analyses presented in this report make use of events recorded by the SD array from 1 
January 2004 to 31 December 2011, with zenith angles less than 55°. To ensure good angle and 
energy reconstructions, each event must satisfy a fiducial cut requiring that the elemental cell of the 
event (that is, the all six neighbours of the water-Cherenkov dete ctor with the highest signal) was 
active when the event was recorded (jAuger Collaboration 2010bl ). Based on this fiducial cut, and 
accounting for unavoidable periods of array instability reducing slightly the duty cycle, the total 
geometric exposure corresponding to the data set considered in this report is 23,520 km 2 yr sr. This 
geometric exposure applies to ene rgies at which the SD array operates with full detection efficiency, 
that is, to energies above 3 EeV (jAuger Collaboration 2010bl ). 



The event direction is determined following the procedure described in Ref. (IBonifazi et al. 20081 ) 
At the lowest energies observed, the angular resolution of the SD is about 2.2°, and reaches ~ 1° 
at the highest energies (jBonifazi et al. 20091 ). This is sufficient to perform searches for large-scale 
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anisotropies. 



The energy estimation of each event is primarily based on the measurement of the signal at a 
reference distance of 1000 m, 5(1000), referred to as the shower size. For a given energy, the shower 
size is a function of the zenith angle due to the rapid increase of the slant depth which induces an 
attenuation of the electromagnetic component of the showers. To account for this attenuation, the 
relationship between the observed 5(1000) and the one that would have been measured had the 
shower a rrived at a z enith angle 38° is derived in an empirical way, using the constant intensity cut 
method (jHersil 196ll ) . To convert 538° i n t° energy, a calibration curve i s used, based on events mea - 



sured simultaneously by the SD array and the fluorescence telescopes ([Auger Collaboration 20081 ). 
since these telescopes indeed provide a calorimetric measurement of the energy. The statistical 
uncertainty of this energy esti mation amounts to about 15 %, while the absolute energy scale has a 
systematic uncertainty of 22% (jAuger Collaboration 2 008). 



3. Control of the event counting rate 

The control of the event counting rate is critical in searches for large scale 
anisotropies. Due to the steepness of the energy spectrum, any mild bias in the estimate of the 
shower energy with time or incident angles can lead to significant distortions of the event counting 
rate. The procedure followed to obtain an unbiased estimate of the shower energy is described 
in this section. This procedure consists in correcting measurements of shower sizes, 5(1000), for 
the influences of weather effects and the geomagnetic field before the conversion to 538° using the 
constant intensity method. Then, the conversion to energy is applied. 



3.1. Influence of atmospheric conditions on shower size 
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Table 1: Coefficients a p , (3 P and a p used to correct shower s izes for atmospheric effects on shower develop- 
ment, in bins of scc6>. From Ref. ( Auger Collaboration 20091) . 



The energy estimator of the showers recorded by the SD array is provided by the signal at 
1000 m from the shower core, 5(1000). For any fixed energy, since the development of extensive 
air showers depends on the atmospheric pressure P and air density p, the corresponding 5(1000) 
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is sensitive to variations in pressure and air density. Systematic variations with time of 5(1000) 
induce variations of the event rate that may distort the real dependence of the cosmic ray intensity 
with right ascension. To cope with this experimental effect, the observed shower size 5(1000), 
measured at the actual density p and pressure P, is related to the one S^ m, (T000) that would have 



been measured at reference values po and Pq (jAuger Collaboration 20091 ) : 



W1000) = [1 - a P (9)(P - P ) - a p (6)(p d - Po ) - P p (9)(p - p d )] 5(1000). (1) 

The reference values are chosen as the average values at Malargiie (i.e. po = 1.06 kg m -3 and 
Po = 862 hPa). pd denotes here the average daily density at the time the event was recorded. The 
measured coefficients a p , (3 P and ap - given in Table Q]- give the influence on the shower sizes of 
the air density (and thus temperature) at long and short time scales on the Moliere radius (and 
hence the lateral profiles of the showers) and of the pressure on the longitudinal development of air 
showers, respectively. 

Applying these corrections to the energy assignments of showers allows us to cancel spurious 
variations of the event rate in right ascension, whose typical amplitudes amount to a few per 
thousand when considering data sets collected over full years. 



3.2. Influence of the geomagnetic field on shower size 

The trajectories of charged particles in extensive air showers are curved in the Earth's magnetic 
field, resulting in a broadening of the spatial distribution of particles in the direction of the Lorentz 
force. As the strength of the geomagnetic field component perpendicular to any arrival direction 
depends on both the zenith and azimuthal angles, the small changes of the density of particles at 
ground induced by the field break the circular symmetry of the lateral spread of the particles and 
thus induce a dependence of the shower size 5(1000) at a fixed energy in terms of the azimuthal 
angle. Due to the steepness of the energy spectrum, such an azimuthal dependence translates into 
azimuthal modulations of the estimated cosmic ray event rate at a given 5(1000). To eliminate 
these effects, the observed shower size 5(1000) is related to the one that woul d have been observed 



in the absence of geomagnetic field 5 9eom (1000) (jAuger Collaboration 2011bl ) : 



5 5eom (1000) = [l - g x cos~ 92 (9) sin 2 (u, b)J 5(1000), (2) 

where g± = (4.2 ± 1) 10~ 3 , <?2 = 2.8 ± 0.3, and u and b = B/||B|| denote the unit vectors in the 
shower direction and the geomagnetic field direction, respectively. At a zenith angle 9 = 55°, the 
amplitude of the asymmetry in azimuth already amounts to ~ 2%, which is why we restrict the 
present analysis to zenith angles smaller than this value. Carrying out these corrections is thus 
critical for performing large scale anisotropy measurements in declination. 
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3.3. Prom shower size to energy 

Once the influence on 5(1000) of weather and geomagnetic effects are accounted for, the 
dependence of 5(1000) on zenith angle due to the attenuation of t he shower and geometrical effects 



is extracted from the data using the constant intensity cut method (jAuger Collaboration 20081 ) . The 
attenuation curve CIC{6) is fitted with a second order polynomial in x = cos (0) - cos 2 (38°) : 
CIC{6) = 1 + ax + bx 2 . The angle 38° is chosen as a reference to convert 5(1000) to 538° = 
5(1000)/CTC(#). 538° m &y be regarded as the signal that would have been expected had the 
shower arrived at 38°. The values of the parameters a = 0.94 ± 0.03 and b = —0.95 ± 0.05 are 
deduced for 538° = 22 VErvfj], that corresponds to an energy of about 4 EeV - just above the 
threshold energy for full efficiency. The differences of these parameters with respect to previous 
reports will be discussed in section [6l 

Finally, the sub-sample of events recorded by both the fluorescence telescopes and the SD 
array is used to establish the relationship between the energy reconstructed with the fluorescence 
telescopes Epn and 538° : Epp, = AS^ a . The resulting parameters from the data fit are A = 
(1.68 ±0.05) x IP -1 Ee V and B = 1.030 ±0.009, in good agreement with the recent report given in 



Ref. (jPesce et al. 201 ll ). The energy scale inferred from this data sample is applied to all showers 



detected by the SD array. 



4. Directional exposure of the Surface Detector array above 1 EeV 

The directional exposure u of the Observatory provides the effective time-integrated collecting 
area for a flux from each direction of the sky H, in units km 2 yr. For energies below 3 EeV, it is 
controlled by the detection efficiency e for triggering. This efficiency depends on the energy E, the 
zenith angle 9, and the azimuth angle tp. Consequently, the directional exposure of the Observatory 
is maximal above 3 EeV, and it is smaller at lower energies where the detection efficiency is less 
than unity. 

In this section we show in a comprehensive way how the directional exposure of the SD array 
is obtained as a function of the energy. We first explain how the slightly non-uniform exposure 
of the sky in sidereal time can be accounted for in the search for anisotropies (section 14. 1 f) . In 
section 14.21 we empirically calculate the detection efficiency as a function of the zenith angle and 
deduce the exposure below the full efficiency energy (3 EeV). In section [4731 we discuss the azimuthal 
dependence of the efficiency due to the geomagnetic effects, introduce the corrections due to the tilt 
of the array in section 14.41 and the corrections due to the spatial extension of the array in section 14.51 



A vertical equivalent muon, or VEM, is the expected signal in a surface detector crossed by a muon traveling 
vertically and centrally to it. 

2 In other contexts such as the determination of the energy spectrum for instance, the term "exposure" refers to 
the total exposure integrated over the celestial sphere, in units km 2 yr sr. 
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and show that the influence of weather effects is negligible on the detection efficiency between 1 and 
3 EeV in section 14.61 Finally we give in section 14.71 some examples of our fully corrected exposure 
at several energies. 



4.1. From local to celestial directional exposure. 

The choice of the fiducial cut to select high quality events allows the precise determination o f 



the geometric directional aperture per cell as a ce n(9) = 1.95 cos 9 km 2 (jAuger Collaboration 2010bl ) 



It also allows us to exploit the regul arity of the array for obtaini ng its geometric directional aper- 



ture as a simple multiple of a ce \\(9) (jAuger Collaboration 2010bl ). The number of elemental cells 



n ceii(t) is accurately monitored every second at the Observatory. To search for celestial large scale 
anisotropics, it is mandatory to account for the modulation imprinted by the variations of n ce n(t) 
in the expected number of events at the sidereal peri odicity T s id- Within each sidereal day, and in 



the same way as in Ref. (jAuger Collaboration 2011al ). we denote by a the local sidereal time and 



express it in hours or in radians, as appropriate. For practical reasons, a is chosen so that it is 
always equal to the right ascension of the zenith at the centre of the array. As a function of a , 
the total number of elemental cells A^ei^a ) and its associated relative variations AA^ ce n(a°) are 
then obtained from : 

N cc ii(a°) = Vn ceU («° + jT sid ), AN ccn (a°) = (3) 



with (iV ce ii) a o = l/T si( i J^ sid da°N ce n(a°). In the same way as in Ref. ( Auger Collaboration 2011al ). 
the small modulation of the expected number of events in right ascension induced by those variations 
will be accounted for by weighting each event k with a factor inversely proportional to AA^ ce n(a^) 
when estimating the anisotropy parameters in section [5j Placing such time dependences in the 
event weights allows us to remove the modulations in time imprinted by the growth of the array 
and the dead times for each detector. 

At any time, the effective directional aperture of the SD array is controlled by the geometric one 
and by the detection efficiency function e(9, (p, E). For each elemental cell, the directional exposure 
in celestial coordinates is then simply obtained through the integration over local sidereal time of 
x^\cP) x a cc \\(9) x e(9,ip,E), where x^ l \a°) is the operational time of the cell (i). Actually, since 
the small modulations in time imprinted in the event counting rate by experimental effects will be 
accounted for by means of the weighting procedure just described when searching for anisotropics, 
the small variations in local sidereal time for each x^ l \a°) can be neglected in calculating oj. The 
zenith and azimuth angles are related to the declination and the right ascension through : 

cos 9 = sin<5 sin£ S ite + cos (5 cos-£ S ite cos (q — a ), 

cos 5 sin 4ite cos (a — a°) — sin 5 cos £ s ite / , \ 

tan^ = \ >— (4) 

cos o sin [a — cr) 
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with ^ s ite the mean latitude of the Observatory. Since both 9 and ip depend only on the difference 
a — a , the integration over a can then be substituted for an integration over the hour angle 
a' = a — q° so that the directional exposure actually does not depend on right ascension when the 
are assumed local sidereal time independent : 

r 2Ah 

da' a cell (9(a', 5)) e(9(a', 6) , <p(a' , 8) , E) . 



"cell 

i=i Jo 



(5) 



Above 3 EeV, this integration can be performed analytically (jSommers 20011 ). Below 3 EeV, the 



non-saturation of the detection efficiency makes the directional exposure lower. The next sections 
are dedicated to the determination of e(9,ip,E). 



4.2. Detection efficiency 

To determine the detection efficiency function, a natural method would be to generate show- 
ers by means of Monte-Carlo simulations and to calculate the ratio of the number of triggered 
events to the total simulated. However, there are discrepancies in the predictions of the hadronic 
intera ction model regar ding the number of muons in shower simulations and what is found in our 



data (jEngel et al. 20071 ). This prevents us from relying on this method for obtaining the detection 



efficiency to the required accuracy. 

We adopt here instead an empirical approach, based on the quasi-invariance of the zenithal 
distribution to large scale anisotropics for zenith angles less than ~ 60° and for any Observatory 
whose latitude is far from the poles of the Earth. For full efficiency, the distribution in zenith angles 
dN/d6 is proportional to sin 9 cos 9 for solid angle and geometry reasons, so that the distribution 
in d-ZV/dsin 2 9 is uniform. Consequently, below full efficiency, any significant deviation from a 
uniform behaviour in the dN/d sin 2 9 distribution provides an empirical measurement of the zenithal 
dependence of the detection efficiency. The quasi-invariance of d N/ d sin 2 9 to large scale anisotropics 
is demonstrated in Appendix A. 

Based on this quasi-invariance, the detection efficiency averaged over the azimuth can be 
estimated from : 

1 dN(sin 2 9,E) 

where the notation (•) stands for the average over tp and the constant Af is the number of events 
that would have been observed at energy E and for any sin 2 9 value in case of full efficiency 
for an energy spectrum dN/dE = 40 (E/ EeV)~ 3 - 27 km" " 2 yr 1 sr 1 EeV 1 - as measured between 



1 and 4 EeV (jAuger Collaboration 2010al ). Consequently, for each zenith angle, this empirical 
measurement of the efficiency provides an estimate relative to the overall spectrum of cosmic rays. 
In particular, since it is applied to all events detected at energy E without distinction based on the 
primary mass of cosmic rays, this technique does not provide the mass dependence of the detection 
efficiency. For that reason, the anisotropy searches reported in section [5] pertain to the whole 
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Fig. 1. — Detection efBciency averaged over the azimuth as a function of sin 2 6 at different energies, 
empirically measured from the data. 



population of cosmic rays, whether this population consists of a single primary mass or a mixture 
of several elements. 

Results are shown in Fig. [T] for four different energies!!. At 4 EeV, a uniform behaviour around 
1 is observed, though quite noisy due to the reduced statistics. This uniform behaviour is consistent 
with full efficiency at this energy, as expected. Note that some values are greater than 1 for energies 
close or higher than 3 EeV, because of the empirical way of measuring the efficiency relative to the 
overall spectrum of cosmic rays. At 2 EeV, a loss of efficiency is observed for vertical showers due 
to the attenuation of the electromagnetic component of the showers. Up to ~ 40°, the detection 
efficiency steadily increases because the projected area of showers at ground gets larger with zenith 
angle. Above ~ 40°, the rapid increase of the slant depth makes then the attenuation of the 
electromagnetic component stronger, but the muonic component of showers becomes dominant and 
ensures a high detection efficiency. At lower energies, the number of muons is, in contrast, too low 
to impact significantly on the detection efficiency above ~ 40° — 45°, so that a clear decrease is 
observed at high zenith angles. In the following, we use parameterisations obtained by fitting each 
distribution with a fourth-order polynomial function in sin 2 6, which is sufficient to reproduce the 
main details as illustrated in Fig. [TJ 



3 To get the detection efficiency at a single energy E, events are actually selected in narrow energy bins around E. 
In addition, to account for the energy spectrum in E~ 3 ' 27 in this energy range, each event is weighted by a factor 

,7,3.27 
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Fig. 2. — Left : Dependence of the detection efficiency on azimuth for = 55° and E = 1 EeV, due 
to geomagnetic effects. Right : Maximal contrast of the azimuthal modulation of the detection efficiency 
induced by geomagnetic effects as a function of the zenith angle. 



4.3. Geomagnetic effects below full efficiency 

In addition to the effects on the energy determination presented in section 13.21 geomagnetic 
effects also affect the detection efficiency for showers with energies below 3 EeV. This is because 
under any incident angles (9,<p), a shower with an energy E triggers the SD array with a probability 
associated with its size which is a function of azimuth because of the geomagnetic effects @ : 
Ex(l + A(0, ip)) B . Above 1 EeV, this effect is in fact the main source of azimuthal dependence of 
the detection efficiency so that to first order in A(0, </>), e(6,<p,E) can be estimated as : 

1 dN(sm 2 0,E(l + A(e,y)) B ) 
<9 ^ E) = N 

* ^E X + B -^ d ^^. (7) 

The correction to the detection efficiency induced by geomagnetic effects, and in particular the 
azimuthal dependence, is thus straightforward to implement from the knowledge of {e(6,<p,E)} . 
An example of such an azimuthal dependence is shown in the left panel of Fig. [21 for E = 1 EeV 
and 9 = 55°. The modulation reflects the one due to the energy determination : the detection 
efficiency is lowered in the directions where the uncorrected energies are under-estimated due to 
geomagnetic effects, and the efficiency is higher where energies are over-estimated. The maximal 
contrast of such azimuthal modulations is displayed in the right panel as a function of the zenith 
angle, for three different energies. At 2 EeV, the amplitude slightly increases up to ~ 35°, staying 
below ~ 0.1%, and then decreases and even cancels due to the saturation of the detection efficiency. 



4 Here, the shorthand notation A(6, if) stands for g\ cos 92 (6) 



sin 2 (u, b) — < sin 2 (u, b) 



The energy E X (1 + 



A(9, ip)) B is actually the one that would have been obtained without correcting for geomagnetic effects. 
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Fig. 3. — Colour-coded altitude (a.s.l.) of the water- Cherenkov detectors. 

In contrast, when going down in energy, the relative amplitude largely increases with the zenith 
angle due to the increase of the derivative term, reaching ~ 1.7% for 9 = 55° and E = 1 EeV. 



4.4. Tilt of the array 

The altitudes above sea level of the water-Cherenkov detectors are displayed in Fig.[3]in colour 
coding. The coordinates are in a Cartesian system whose origin is defined at the "centre" of the 
Observatory site. The Andes ridge building up in the western and north-western direction can be 
seen. A slightly tilted SD array gives rise to a small azimuthal asymmetry, and consequently slightly 
modifies the directional exposure with respect to Eqn. [5] through small changes of the geometric 
directional aperture. This modification is twofold : the tilt changes the geometric factor (cos#) of 
the projected surface under incidence angles (9,<p); and also induces a compensating effect below 
full efficiency by slightly varying the detection efficiency with the azimuth angle <p. 

Denoting the normal vector to each elemental cell, the geometric directional aperture per 
cell is not any longer simply given by cos 9 but now depends on both 9 and ip : 

a^ u (9,ip) = 1.95 n-n^ ~ 1.95 [1 + C W tan(9cos (<p - iff)] cos6>, (8) 

where and ip^ are the zenith and azimuth angles of n^. It is actually this latter expression 
a ce \\ which has to be inserted into Eqn.[5]to calculate the directional exposure. Overall, the average 
tilt of the SD array is £ eff ~ 0.2°, and induces a dipolar asymmetry in azimuth with a maximum 
in the downhill direction (pf^ ~ 0° and with an amplitude increasing with the zenith angle as 
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~ 0.3% tan 6. 

Below 3 EeV, the tilt of the array induces an additional variation of the detection efficiency 
with azimuth. This is because the effective separation between detectors for a given zenith angle 
depends now on the azimuth. Since, for a given zenith angle, the SD array seen by showers coming 
from the uphill direction is denser than that for those coming from the downhill direction, the 
detection efficiency is higher in the uphill direction. Parameterising the energy dependence of e as 
E 3 /(E 3 +Eq 5 ), we show in Appendix B that the change in the detection efficiency can be estimated 
as : 

E 3 (El 5 -E$f(9,y)) 
(EZ + El 5 )(E3 + E$f(6,<p)) 

where Eq}£(6,(p) is related to £0.5 through : 



E$%(0,<p) * £ .5 x [1 + C eff tan cos (if - ^f)] 3/2 . (10) 

Around 1 EeV, this correction tends to compensate the pure geometrical effect described above, 
and even overcompensates it at lower energies. 



4.5. Spatial extension of the array 

This spatial extension of the SD array is such that the range of latitudes covered by all cells 
reaches ~ 0.5°. This induces a slightly different directional exposure between the cells located at 
the northern part of the array and the ones located at the southern part. This spatial extension can 
be accounted for to calculate the overall directional exposure using the cell latitudes instead 
of the mean site one in the transformations from local to celestial angles in Eqn. [U 



4.6. Weather effects below full efficiency 

In the same way as geomagnetic effects, weather effects can also affect the detection efficiency 

for showers with energies below 3 EeV. However, above 1 EeV, we have shown in (jAuger Collaboration 2011al ) 
that as long as the analysis covers an integer number of years with almost equal exposure in every 
season, the amplitude of the spurious modulation in right ascension induced by this effect is small 
enough to be neglected when performing anisotropy analyses at the present level of sensitivity. 
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Fig. 4. — Directional exposure lo(5,E) as a function of the declination 5, for three different energies. 



4.7. Final estimation of the directional exposure - Examples at some energies 

Accounting for all effects, the final expression to calculate the directional exposure is slightly 
modified with respect to Eqn. [5] : 

"cell „24h 

u(6,E) = Y;x® do , a<2i(0,yO [e(e,<p,E) + Ae m (0,<P,E)], (11) 

where both 6 and tp depend on a', 5 and The resulting dependence on declination is displayed 
in Fig. 2] for three different energies. Down to 1 EeV, the detection efficiency at high zenith angles 
is high enough that the equatorial south pole is visible at any time and hence constitutes the 
direction of maximum of exposure. For a wide range of declinations between ~ —89° and ~ —20°, 
the directional exposure is ~ 2, 500 km 2 yr at 1 EeV, and ~ 3, 500 km 2 yr for any energy above full 
efficiency. Then, at higher declinations, it smoothly falls to zero, with no exposure above ~ 20° 
declination. 

The average expected number of events within any solid angle and any energy range can be 
recovered by integrating the directional exposure over the solid angle considered and the cosmic ray 
energy spectrum in the corresponding energy range. Note that the rapid variation of the exposure 
close to the South pole on an angular scale of the order of the angular resolution has no influence on 
the event counting rate, due to the quasi-zero solid angle in that particular direction. Consequently, 
though the exposure around the South pole could be affected by small changes of the detection 
efficiency around 9 = 55°, the results presented in next sections are on the other hand not affected 
by the exact value of the exposure for declinations a few degrees away from the South pole. 
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5. Searches for large scale patterns 

5.1. Estimates of spherical harmonic coefficients 

Any angular distribution over the sphere <£(n) can be decomposed in terms of a multipolar 
expansion : 

i 

^ = a emYe m (n), (12) 

£>0 m =-l 

where n denotes a unit vector taken in equatorial coordinates. The customary recipe to extract 
each multipolar coefficient makes use of the completeness relation of spherical harmonics : 



a im = / cM $(n)y €m (n), (13) 

J Am 

where the integration is over the entire sphere of directions n. Any anisotropy fingerprint is encoded 
in the a^ m spherical harmonic coefficients. Variations on an angular scale of radians contribute 
amplitude in the I ~ 1/0 modes. 

However, in case of partial sky coverage, the solid angle in the sky where the exposure is 
zero makes it impossible to estimate the multipolar coefficients a^ m in this way. This is because 
the unseen soli d angle prevent s one from making use of the completeness relation of the spher- 



ical harmonics (jSommers 200 ll ) . Since the observed arrival direction distribution is in this case 
the combination of the angular distribution <3?(n) and of the directional exposure function w(n), 
the integration performed in Eqn. [13] does not allo w any longer the extrac tion of the multipolar 



coefficients of <3?(n), but only the ones of w(n) 3>(n) (|Billoir Delignv 20081 ) rl: 



Hm 



dn w(n)$(n)r £m (n) 
Ail 



£'>0m' 



a Vm' I dfiw(n)F f - m -(n)F fm (n). (14) 
,_ „, J Ail 



Formally, the d£ m coefficients appear related to the be m ones through a convolution such that 
bim = Yle'>o Em '=-i'\^tim a t'm' ■ The matrix K, which imprints the interferences between modes 
induced by the non-uniform and partial coverage of the sky, is entirely determined by the directional 
exposure. The relationship established in Eqn. Q3]is valid for any exposure function u;(n). 

Meanwhile, the observed arrival direction distribution, dN(n)/dQ, provides a direct estimation 
of the b£ m coefficients through (hereafter, we use an over-line to indicate the estimator of any 
quantity) : 

b tm = [ dn Y em (n), (15) 

J An asz 



5 To cope with the unseen solid angl e, another approach ma kes use of orthogonal functions of increasing multipo- 
larity, tailored to the exposure u itself ( Billoir fc Delignv 2 008). This method would yield similar accuracies. 
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where the distribution diV(n)/df2 of any set of ./V arrival directions {ni,...,njv} can be modelled 
as a sum of Dirac functions on the sphere. Then, if the multipolar expansion of the angular 
distribution 3>(n) is bounded to £ max , that is, if the <&(n) has no higher moments than ^ maX ; the 
first b£ m coefficients with t < £ max are related to the non- vanishing ai m by the square matrix -?Q max 
truncated to £ max - Inverting this truncated matrix allows us to recover the underlying ag m from 
the measured bi m (with i < £ max ) : 

a im = Yl E [^U^'W- ( 16 ) 
e'=om'=-e> 

In the case of small anisotropies (\ae m \/aoo C 1), the resolution on each recovered ag m coefficient 

0.5 



is proportional to (AK^J$£\ (jBilloir & Delignv 200^ : 



([^Jt^oo) ' 5 - (17) 



The dependence on iraax of the coefficients of induces an intrinsic indeterminacy of each 

recovered coefficient ag m as £ max is increasing. This is nothing else but the mathematical translation 
of it being impossible to know the angular distribution of cosmic rays in the uncovered region of 
the sky. 

Henceforth, we adapt this general formalism to the search for anisotropies in Auger data in 
different energy intervals. We assume that the energy dependence of the angular distribution 
of cosmic rays is smooth enough that the multipolar coefficients can be considered constant for 
any energy E within a narrow interval AE. The directional exposure is hereafter considered as 
independent of the right-ascension, as defined in section [H Within an energy interval AE, the 
expected arrival direction distribution thus reads : 

^^occD(5)^^ a em Y em (n), (18) 

e>om=-e 

where oj(S) is the effective directional exposure for the energy interval AE. For convenience, this 
latter function is normalised such that : 
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f dEE~^oj(5,E) 




AE 


max 


[ dEE-^oj(S,E) 
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Jae 



(19) 



with 7 the spectral index in the considered energy range. This dimensionless function provides, 
for any direction on the sky, the effective directional exposure in the energy range AE at that 
direction, relative to the largest directional exposure on the sky. This is actually the relevant 
quantity which enters into Eqn. Q3] for the analyses presented below. Note that for a directional 
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exposure independent of the right ascension, the coefficients are proportional to 5™ - i.e. 

different values of m are not mixed in the matrix. The observed arrival direction distribution, 
diV(n)/dr2, is here modelled as a sum of Dirac functions on the sphere weighted by the factor 
/S.N~2\{c*k) f° r each event recorded at local sidereal time as described in section H~T1 to correct 
for the slightly non-uniform directional exposure in right ascension. In this way, the integration in 
Eqn. [T4l yields to : 

£j AiV cc ii(a£) 

The multipolar coefficients at m are then recovered by means of Eqn. [TBJ Given the exposure 
functions described in section 0J the resolution on each recovered coefficient, encoded in Eqn. 1171 is 
degraded by a factor larger than 2 each time £ max is incremented by 1. This prevents the recovery 
of each coefficient with good accuracy as soon as l max ^ 3, since, for ^ max — 3 for instance, our 
current statistics would only allow us to probe dipole amplitudes at the 10% level. Consequently, 
in the following, we restrict ourselves to reporting results on individual coefficients obtained when 
assuming a dipolar distribution (£ max = 1) and a quadrupolar distribution (^ max = 2). Meanwhile, 
due to the interferences between modes induced by the non-uniform and partial sky coverage, it 
is important to stress again that each multipolar coefficient recovered under the assumption of a 
particular bound ^ max might be biased if the underlying angular distribution of cosmic rays is not 
bounded to £ max . Given the directional exposure functions considered in this study, this effect can 
be important only if the angular distribution has in fact significant moments of order £ max + 1. 



5.2. Searches for dipolar patterns 

As outlined in the introduction, a measurable dipole is regarded as a likely possibility in many 
scenarios for the origin of cosmic rays at EeV energies. Assuming that the angular distribution of 
cosmic rays is modulated by a pure dipole, the intensity $(n) can be parameterised in any direction 
n as : 

$(n) = — ( 1 + r d-n], (21) 



4-7T 

where d denotes the dipole unit vector. The dipole pattern is here fully characterised by a dec- 
lination 5d, a right ascension ct^, and an amplitude r corresponding to the maximal anisotropy 
contrast : 

^max ™min / nn \ 

r = ^ ^ • I 22 ) 

<J> _l_ <T> v ' 

Y max * ^min 

The estimation of these three coefficients is straightforward from the estimated spherical har- 
monic coefficients a\ m : r = [3(af + a\ 1 + a^_ 1 )] ' 5 /aoo> <5 = arcsin (\/3aio/ooo^)) an d a = 
arctan (ai-i/an). Uncertainties on r, 5 and a are obtained from the propagation of uncertain- 
ties on each recovered a\ m coefficient (cf Eqn. fT7|) . Under an underlying isotropic distribution, 
and for an axisymmetric directional exposure around the axis defined by the North and South 
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5=20° 




Fig. 5. — Left : Reconstructed amplitude of the dipole as a function of energy. The dotted line stands for the 
99% C.L. upper bounds on the amplitudes that would result from fluctuations of an isotropic distribution. 
Right : Reconstructed declination and right-ascension of the dipole with corresponding uncertainties, as a 
function of energy, in azimuthal projection. 



equatorial poles, the probability density function of r is given by ([Auger Collaboration 201 lbl ) : 

(23) 



PR\r) 



erfi 



^Ja\ — a 2 r 



GO, 



V2 



exp 



where ern(z) = eif(iz)/i, a = V3au/aoo, and a z = \/3cio/aoo- The probability Pr(> r) that an 
amplitude equal or larger than r arises from a statistical fluctuation of an isotropic distribution is 
then obtained by integrating pn above r : 



Pr{> f) = erfc 



VV2, 



+ erfi 



y 1 o 2 — a 2 r 
crcr z ^2 



cxp 



2a 2 



(24) 



The reconstructed amplitudes r{E) and corresponding directions are shown in Fig. [5] with 
the associated uncertainties, as a function of the energy. The directions are drawn in azimuthal 
projection, with the equatorial South pole located at the centre and the right-ascension going from 
to 360° clockwise. In the left panel, the 99% C.L. upper bounds on the amplitudes that would result 
from fluctuations of an isotropic distribution are indicated by the dotted line (i.e. the amplitudes 
Tgg(E) such that Pr(> rgg(E)) = 0.01). One can see that within the statistical uncertainties, there 
is no strong evidence of any significant signal. 

The reconstructed declinations 5 and right ascensions a are shown separately in Fig [6l Both 
quantities are expected to be randomly distributed in case of independent samples whose par- 
ent distribution is isotropic. In our previous report on first harmonic analysis in right ascen- 



sion (jAuger Collaboration 201 lal ) , we pointed out the intriguing smooth alignment of the phases in 
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Fig. 6. — Reconstructed decli nation (left) and right ascen sion (right) of the dipole as a function of energy. 
The smooth fit to the data of (| Auger Collaboration 201 lal ) is shown as the dashed line in the right panel : a 
consistent smooth behaviour is observed using the analysis presented here and applied to a data set containing 
two additional years of data. 
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Fig. 7. — Significance sky maps in four independent energy bins. The maps are smoothed using an angular 
window with radius = 1 radian, to exhibit any dipolar-like structures. The directions of the reconstructed 
dipoles are shown with the associated uncertainties. The galactic plane and galactic center are also depicted 
as the dotted line and the star. 
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Fig. 8. — Left : Amplitude of the dipole for two energy intervals : 1 < i?/[EeV] < 4 and E > 4 EeV. Right : 
Amplitude of the dipole as a function of energy thresholds. The dotted lines stand for the 99% C.L. upper 
bounds on the amplitudes that could result from fluctuations of an isotropic distribution. 



right ascension as a function of the energy, and noted that such a consistency of phases in adjacent 
energy intervals is expected to manifest with smaller number of events than those required for 
the detection of amplitudes standing-out significantly above the background noise in case of a real 
underlying anisotropy. This motivated us to design a prescription aimed at establishing at 99% 
C.L. w hether this consistency in ph ases is real, using the exact same analysis as the one reported 
in Ref. (1 Auger Collaboration 2011al ), The prescribed test will end onc e the total exposure since 25 
June 2011 is 21,000 km 2 yr sr. The smooth fit to the data of Ref. (jAuger Collaboration 201 lal ) 
is shown as a dashed line in the right panel of Fig [6j restricted to the energy range considered 
here. Though the phase between 4 and 8 EeV is poorly determined due to the corresponding di- 
rection in declination pointing close to the equatorial south pole, it is noteworthy that a consistent 
smooth behaviour is observed using the analysis presented here and applied to a data set containing 
two additional years of data. It is also interesting to see in the left panel that all reconstructed 
declinations are in the equatorial southern hemisphere. 

For completeness, significance sky maps are displayed in Fig. [7] in equatorial coordinates and 
using a Mollweide projection, for the four energy ranges. The galactic plane and galactic center 
are also de picted as the d otted line and the star. Significances are calculated using the Li and Ma 
estimator (|Li &: Ma 19831 ). This widely used estimator of significance, S, properly accounts for the 
fluctuations of the background and of an eventual signal in any angular region searched |§. If no 
signal is present, the variable S is nearly normally distributed even for small count numbers, so that 
positive values of S can be interpreted as the number of standard deviations of any excess in the 



6 The parameter olm in the expression of the Li & Ma significance, expressing the expected ratio of the count 
numbers between the angular region searched (the on-region) and any background region if there is no signal in the 
on-region, is here taken as the ratio between the expected number of events in the on-region and the total number of 
events in the energy range considered. 
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sky. As well, for negative values of S, —S can be interpreted as the number of standard deviations 
of any deficit in the sky. The maps show the overdensities obtained in circular windows of radius 
= 1 radian, to better exhibit possible dipolar-like structures. The directions of the reconstructed 
dipoles are also shown, with their associated uncertainties (thick circles). 

Finally, since some consistency is observed both in declination and right ascension as a func- 
tion of energy, the use of larger energy intervals and/or energy thresholds may help to pick up a 
significant signal above the background level. The amplitudes of the dipole are shown in Fig. [8] for 
two energy intervals (1 < £7/ [EeV] < 4 and E > 4 EeV) and as a function of energy thresholds. 
This does not provide any further evidence for significant anisotropics. 



5.3. Searches for quadrupolar patterns 

Any excesses along a plane would show up as a prominent quadrupole moment. Such excesses 
are plausible for instance at EeV energies in case of an emission of light EeV-cosmic rays from 
sources preferentially located in the galactic disk, or at higher energies from sources preferentially 
located in the super-galactic plane. Consequently, a measurable quadrupole may be regarded as an 
interesting outcome of an anisotropy search at ultra high energies. 

Assuming now that the angular distribution of cosmic rays is modulated by a dipole and a 
quadrupole, the intensity <3?(n) can be parameterised in any direction n as : 

cD(n) = ^(l + rd.n + ^ Qi^n^j , (25) 

i,j 

where Q is a traceless and symmetric second order tensor. Its five independent components are 
determined in a straightforward way from the £ = 2 spherical harmonic coefficients <Z2 m . Denoting 
by A + ,Ao,A_ the three eigenvalues of Q/2 (A+ being the highest one and A_ the lowest one) and 
q+, qo, q~ the three corresponding unit eigenvectors, the intensity can be parameterised in a more 
intuitive way as : 



$(n) = 4^ V 1 + r d ' " + A+(q+ ' n) + A ° (q ° ' n) + A ~ (q - ' n) ) ■ (26) 
It is then convenient to define the quadrupole amplitude (3 as : 

^^TT' < 27 > 

In case of a pure quadrupolar distribution (i.e. in the absence of dipole), j3 is nothing else but the 
customary measure of maximal anisotropy contrast : 

r = o p = A+ ~ A - = $max ~^ min . (28) 
2 + A+ + A_ 

min 
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Fig. 9. — Amplitudes of the dipolar (top) and quadrupolar moments (middle and bottom) as a function 
of energy using a multipolar reconstruction up to Cax = 2, for two different binnings (left and right). In 
each panel, the dotted lines stand for the 99% CL. upper bounds on the amplitudes that could result from 
fluctuations of an isotropic distribution. 



Hence, any quadrupolar pattern can be fully described by two amplitudes (/3, A+) and three angles : 
(<5+,a+) which define the orientation of q+ and («-) which defines the direction of q_ in the 
orthogonal plane to q + . The third eigenvector qo is orthogonal to q + and q_, and its corresponding 
eigenvalue Aq is such that the traceless condition is satisfied : A + + A„ + Aq = 0. Though the 
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probability density functions of the estimated quadrupole amplitudes A+) can be in principle 
calculated in the same way as in the case of the estimated dipole amplitude (r), expressions are 
much more complicated to obtain even semi-analytically and we defer hereafter to Monte-Carlo 
simulations to tabulate the distributions. 

The amplitudes r(E), \+(E) and /3(E) are shown in Fig. as functions of energy. Dipole 
amplitudes are compatible with expectations from isotropy. Compared to the results on the dipole 
obtained in previous section for ^ max = 1> the sensitivity is now degraded by a factor larger than 2 
as expected from the dependence of the resolution ai m on ^ max (cf Eqn. [T7|) . In the same way as for 
dipole amplitudes, the 99% C.L. upper bounds on the quadrupole amplitudes that could result from 
fluctuations of an isotropic distribution are indicated by the dashed lines. They correspond to the 
amplitudes \ +> gg(E) and f3 99 (E) such that the probabilities P\ + (> X +j gg(E)) and Pb(> figg{E)) 
arising from statistical fluctuations of isotropy are equal to 0.01. Here, both distributions P\ + and 
Pb are sampled from Monte-Carlo simulations. Throughout the energy scan, there is no evidence 
for anisotropy. The largest deviation from isotropic expectations occurs between 2 and 4 EeV, 
where both amplitudes A + and f3 lie just above A + gg and f3 99 . 

6. Additional cross-checks against experimental effects 

6.1. More on the influence of shower size corrections for geomagnetic effects 

Understanding the influence of the shower size corrections for geomagnetic effects is critical to 
get unbiased estimates of anisotropy parameters. Without accounting for these effects, an increase 
of the event rate would be observed close to the equatorial South pole with respect to expectations 
for isotropy, while a decrease would be observed close to the edge of the directional exposure in 
the equatorial Northern hemisphere. This would result in the observation of a fake dipole. A 
convenient way to exhibit this effect is to separate the dipole in two components : the component 
of the dipole in the equatorial plane r±, and the component along the Earth rotation axis, rii. 
While rj_ is expected to be affected only by time-dependent effects, rn is on the other hand the 
relevant quantity sensitive to time-independent effects such as the geomagnetic one. 



AE [EeV] 


r7 corr [%] 


r±[%] 


jflincorr 


f,|[%] 


1-4 


0.9 ±0.3 


0.9 ±0.3 


-2.2 ±0.4 


-1.0 ±0.4 


> 4 


1.8 ±1.0 


2.1 ± 1.0 


-4.1 ± 1.7 


-3.0 ± 1.7 



Table 2: Influence of shower size corrections for geomagnetic effects on the component of the dipole in the 
equatorial plane and on the one along the Earth rotation axis. 

Estimations of r± and rii obtained by accounting or not for geomagnetic effects are given in 
Table El in two different energy ranges. These estimations are obtained from the recovered ~a\ m 
coefficients : rii = v^aio/<%b and r± = + «i-i]°' 5 /ooo- It can be seen that the main effect of 
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the geomagnetic corrections is a shift in rii of about 1.2%. In the energy range 1 < E/\EeV] < 4, 
this shift is significant, fu changing from -2.2% to -1.0% with an uncertainty amounting to 0.4%. 
Above 4 EeV, the net correction is of the same order, though the statistical uncertainties are larger. 
In contrast, r± remains unchanged in both cases, as expected. 

6.2. Eventual energy dependence of the attenuation curve 

In this section, we study to which extent the procedure used to obtain the attenuation curve 
in section 13.31 might influence the determination of the anisotropy parameters. 

To convert the shower size into energy, we explained and applied in section 13.31 the constant 
intensity cut method for showers with 538° > 22 VEM, that is, just above the threshold energy 
for full efficiency. The value of the parameter a obtained in these conditions is consistent within 
the statistical uncertainties with the one previously reported when applying the same constant 
intensity cut method for showers with £38° > 47 VEM. Opposite to this, the value obtained for the 
coefficient b differs by more than 3 standard deviations. Such a difference might be expected from 
both the evolution of the maximum of the showers and from an eventual change in composition 
with energy, but it may also be due to energy and angle-dependent resolutions effects mimicking a 
real evolution with energy. 

With a different attenuation curve, some events would be reconstructed in the adjacent energy 
intervals in an extent which depends on the change of the attenuation curve with zenith angle. For 
that reason, the determination of anisotropy parameters might be altered by this effect. 

Disentangling real evolution of the attenuation curve with energy from resolution effects is 
out of the scope of this paper and will be addressed elsewhere. Here, we restrict ourselves to 
probe the effect that a real energy dependence would have on the determination of anisotropy 
parameters. To do so, we choose to fit the values of the coefficient b obtained for S^s° = 22 VEM 
and £38° = 47 VEM through a linear dependence with the logarithm of 5*38° • Below and above 
these values, the behaviour of b(E) is obtained by extrapolating this energy dependence. In this 
way, the changes in the anisotropy parameters are probed in extreme conditions. 

Repeating the whole chain of analysis with this new attenuation curve, it turns out that the 
reconstructed dipole parameters are only marginally affected by this change, as illustrated in the 
top and middle panels of Fig. [lOl Meanwhile, both reconstructed quadrupole amplitudes in the 
energy interval 2 < E/~EeV < 4 are reduced in such a way that they lie now just below the 99% 
upper bounds for isotropy. Conversely, the amplitudes in the energy interval 1 < E/~EeV < 2 
are slightly increased. Below 4 EeV, the determination of the attenuation curve thus appears to 
bring some systematic uncertainties for determining the quadrupole amplitudes. The two extreme 
extrapolations performed in this analysis (i.e. b constant with the energy or linearly dependent 
with the logarithm of the energy) allows us to bracket the possible values. 
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Fig. 10. — Impact of different sources of systematic uncertainties on the dipole amplitudes (top) and the 
dipole directions and phases (middle) obtained under the assumption £ max = 1, and quadrupole amplitudes 
(bottom) obtained with ^ max = 2, as a function of the energy. The blue bands correspond to the results 
presented in Fig. [5] and Fig. GO 
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6.3. Systematic uncertainties associated to corrections for weather and geomagnetic 

effects 

In section [31 we presented the procedure adopted to account for the changes in shower size due 
to weather and geomagnetic effects. Since the coefficients ap, a p and f3 p in Eqn. [T]were extracted 
from real data, they suffer from statistical uncertainties which may impact in a systematic way the 
corrections made on S'(IOOO), and consequently may also impact the anisotropy parameters derived 
from the data set. Besides, the determination of g\ and 52 in Eqn. [2] is based on the simulation 
of showers. Both the systematic uncertainties associated to the different interaction models and 
primary masses and the statistical uncertainties related to the procedure used to extract g± and 52 
constitute a source of systematic uncertainties on the anisotropy parameters. 

To quantify these systematic uncertainties, we repeated the whole chain of analysis on a large 
number of modified data sets. Each modified data set is built by sampling randomly the coefficients 
ap, a p and j3 p (or g\ and g% when dealing with geomagnetic effects) according to the corresponding 
uncertainties and correlations between parameters through the use of a Gaussian probability distri- 
bution function. For each new set of correction coefficients, new sets of anisotropy parameters are 
then obtained. The RMS of each resulting distribution for each anisotropy parameter is the system- 
atic uncertainty that we assign. Results are shown in Fig. [TUl in terms of the dipole and quadrupole 
amplitudes as a function of the energy. Balanced against the statistical uncertainties in the original 
analysis (shown by the bands), it is apparent that both sources of systematic uncertainties have a 
negligible impact on each reconstructed anisotropy amplitude. 



7. Upper limits and discussion 

From the analyses reported in section 02 upper limits on dipole and quadrupole amplitudes can 
be derived at 99% C.L. (see appendices C and D). All relevant results are summarised in Table [3] 
and Table HI The upper limits are also shown in Fig. II II accounting for the systematic uncertainties 
discussed in the previous section : in the two last energy bins, the upper limits are quite insensitive 
to the systematic uncertainties because all amplitudes lie well within the background noise. 

We illustrate below the astrophysical interest of these upper limits by calculating the anisotropy 
amplitudes expected in a toy scenario in which sources of EeV-cosmic rays are stationary, densely 
and uniformly distributed in the galactic disk, and emit particles in all directions. 

Both the strength and the structure of the magnetic field in the Galaxy, known only approx- 
imately, play a crucial role in the propagation of cosmic rays. The field is thought to contain a 
large scale regular compo nent and a small scale turbulent one, both having a local strength of a 



few microgauss (see e.g. (jBeck 20011 )). While the turbulent component dominates in strength by 
a factor of a few, the regular component imprints dominant drift motions as soon as the Larmor 
radius of cosmic rays is larger than the maximal scale of the turbulences (thought to be in the range 
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AE [EeV] 


N 


r [%} 


5[°] 


op] 


UL [%] 


1 - 2 


360132 


1.0 ±0.4 


-15 ±32 


342 ± 20 


1.5 


2-4 


88042 


1.6 ±0.8 


-46 ± 28 


35 ±30 


2.8 


4-8 


19794 


2.7 ±2.0 


-69 ± 30 


25 ±74 


5.8 


> 8 


8364 


7.5 ±2.5 


-37 ±21 


96 ± 18 


11.4 



Table 3: Summary of the dipolar analysis (£ max = 1) reported in section \5?2\ together with the derived 99% 
C.L. upper limits (UL) on the amplitudes. 



AE [EeV] 


A+ [%] 


n%\ 


UL (A+) [%] 


UL 09) [%] 


1 - 2 


2.0 ±0.7 


1.7 ±0.6 


3.0 


2.9 


2-4 


5.0 ± 1.7 


4.2 ±1.3 


6.3 


6.1 


4-8 


1.6 ±2.0 


1.9 ± 1.8 


10.0 


9.4 


> 8 


4.0 ±3.4 


3.9 ±2.7 


14.5 


13.8 



Table 4: Summary of the quadrupolar analysis (£ max = 2) reported in section [531 together with the derived 
99% C.L. upper limits (UL) on the amplitudes. 



10-100 pc). We adopt in the following a recent parameterisation of the regular component obtained 
by fitting model field geometrie s to Faraday rotation measures of extragalactic radio sources and 



polarised synchrotron emission (jPshirkov et al. 201 ll ). It consists in two different components : a 



disk field and a halo field. The disk field is symmetric with respect to the galactic plane, and is 
described by the widely-used logarithmic spiral model with reversal direction of the field in two dif- 
ferent arms (the so-called BSS-model). The halo field is anti-symmetric with respect to the galactic 



plane and purely toroidal. The detailed parameterisation is given in Ref. (jPshirkov et al. 20111 ) 
(with the set of parameters reported in Table 3). In addition to the regular component, a turbu- 
lent field is generated according to a Kolmogorov power spectrum and is pre-computed on a three 
dimensional grid periodically repeated in space. The size of the grid is taken as 100 pc, so as the 
maximal scale of turbulences, and the strength of the turbulent component is taken as three times 
the strength of the regular one. 

To describe the propagation of cosmic rays with energies E > 1 EeV in such a magnetic field, 
the direct integration of trajectories is the most appropriate tool. Performing the forward tracking 
of particles from galactic sources and recording those particles which cross the Earth is however not 
feasible within a reasonable computing time. So, to obtain the anisotropy of cosmic rays emitted 
from sources uniformly distributed in a disk with a radius of 20 kp c from the galactic centre and with 



a height of ± 100 pc, we adopt a method first proposed in Ref. (jThielheim Langhoff 19681 ) and 
then widely used in the literature. It consists in back tracking anti-particles with random directions 
from the Earth to outside the Galaxy. Each test particle probes the total luminosity along the path 
of propagation from each direction as seen from the Earth. For stationary sources emitting cosmic 
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Fig. 11. — 99% C.L. upper limits on dipole and quadrupole amplitudes as a function of the energy. Some 
generic anisotropy expectations from stationary galactic sources distributed in the disk are also shown, for 
various assumptions on the cosmic ray composition. The fluctuations of the amplitudes due to the stochastic 
nature of the turbulent component of the magnetic held are sampled from different simulation data sets and 
are shown by the bands (see text). 



rays in all directions, the flux expected in a given sampled direction is then proportional to the 
time spent in the source region by the test particles arriving from that direction. 

The amplitudes of anisotropy obviously depend on the rigidity E/Z of the cosmic rays, with Z 
the electric charge of the particles. Since we only aim at illustrating the upper limits, we consider 
two extreme single primaries : protons and iron nuclei. In the energy range 1 < E/~EeV < 20, it is 
unlikely that our measurements on the average position in the atmosp here of the shower maximum 
and the corresponding RMS can be reproduced with a single primary (jAuger Collaboration 2010q ) . 
As well, in the scenario explored here and for a single primary, the energy spectrum is expected 
to reveal a hardening in this energy range, whose origin is different from th e one expecte d if the 
ankle marks the cross-over between galactic and extragalactic cosmic rays (jLinslev 19631 ) or if it 
marks the distortion of a proton-dominated extragalactic spectrum due to e + / e ~ pair production 
of protons with the photons of the cosmic mic rowave background (IHillas 19671 : iBlumenthal 19701 : 
Berezinsky et al. 20061 : iBerezinsky et al. 20041 ). For a given configuration of the magnetic field, 
the exact energy at which this hardening occurs depends on the electric charge of the cosmic rays. 
This is because the average time spent in the source region first decreases as ~ E^ 1 and then 
tends to the constant free escape time as a consequence of the direct escape from the Galaxy. The 
hardening with A7 ~ 0.6 observed at 4 EeV in our measurements of the energy spectrum is not 
compatible with the one expected in this scenario (A7 ~ 1). Nevertheless, the calculation of dipole 
and quadrupole amplitudes for single primaries is useful to probe the allowed contribution of each 
primary as a function of the energy. 

The dipole r and quadrupole A+ amplitudes obtained for several energy values covering the 
range 1 < E/~EeV < 20 are shown in Fig. [TTl To probe unambiguously amplitudes down to the 
percent level, it is necessary to generate simulated event sets with ~ 5 10 5 test particles. Such 
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a number of simulated events allows us to shrink statistical uncertainties on amplitudes at the 
0.5% level. Meanwhile, there is an intrinsic variance in the model for each anisotropy parameter 
due to the stochastic nature of the turbulent component of the magnetic field. This variance is 
estimated through the simulation of 20 sets of 5 10 5 test particles, where the configuration of the 
turbulent component is frozen in each set. The RMS of the amplitudes sampled in this way is 
shown by the bands in Fig. [TTJ While the dipole amplitude steadily increases for iron nuclei, this 
is not the case any longer for protons around the ankle energy. This is because we explore a source 
region uniformly distributed in the disk. Consequently, the image of the galactic plane appears less 
distorted by the magnetic field with increasing energy. This gives rise to an important quadrupolar 
moment which actually turns out to be the main feature of the anisotropy at large scale 0. 

The dipole and quadrupole A + amplitudes obtained here depend on the mo del used to describe 



the galactic magnetic field. We note that recently, a new model was given in Ref. (IFarrar &: Jansson 20121 ). 
providing improved fits to Faraday rotation measures of extragalactic radio sources and polarised 
synchrotron emission observations. However, we tested at a few energies that the results obtained 
are qual itatively in agreemen t with the ones presented in Fig. [TTJ Similar conclusions were given 



in Ref. (jGiacinti et al. 2012Q . where more systematic studies can be found in terms of the field 



strength and geometry. 

Around 1 EeV, there are indications that the cosmic ray composition inclu des a significant light 



component from var i ous measuremen ts of the depth of shower maximum X max (lAuger Collaboration 2010d : 



Abbasi et al. 2010bl : lJui et al. 20111 ). It is apparent that amplitudes derived for protons largely 
stand above the allowed limits. Consequently, unless the strength of the magnetic field is much 
higher than in the picture used here, the upper limits derived in this analysis exclude that the light 
component of cosmic rays comes from galactic stationary sources densely distributed in the galactic 
disk and emitting in all directions. This is in agreement with the absence of any detectable point-like 
sources above 1 EeV that would be indicative of a flux of neutrons produced by EeV-protons th rough 



mainly pion-producing interactions in the source environments ([Auger Collaboration 2012bl ). On 
the other hand, if the cosmic ray composition around 1 EeV results from a mixture containing a 
large fraction of iron nuclei of galactic origin, upper limits can still be respected, or alternatively 
a light component of extragalactic origin would be allowed. Future measurements of composition 
below 1 EeV will come from the low energy extension HEAT now available at the Pierre Auger 



Observatory (jMathes et al. 201 ll ) . Combining these measurements with large scale anisotropy ones 



will then allow us to further understand the origin of cosmic rays at energies less than 4 EeV. 



7 This feature would remain in the case of a radial distribution of sources following the matter in the Galaxy, 
though the dipole amplitude would steadily increase above the ankle energy. 
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8. Summary 

For the first time, a thorough search for large scale anisotropies as a function of both the 
declination and the right ascension in the distribution of arrival directions of cosmic rays detected 
above 1 EeV at the Pierre Auger Observatory has been presented. With respect to the traditional 
search in right ascension only, this search requires the control of additional systematic effects 
affecting both the exposure of the sky and the counting rate of events in local angles. All these 
effects were carefully accounted for and presented in sections [3] and 01 No significant deviation from 
isotropy is revealed within the systematic uncertainties, although the consistency in the dipole 
phases may be indicative of a genuine signal whose amplitude is at the level of the statistical noise. 
The sensitivity accumulated so far to dipole and quadrupole amplitudes allows us to challenge an 
origin of cosmic rays from stationary galactic sources densely distributed in the galactic disk and 
emitting predominantly light particles in all directions. 

Future work will profit from both the increase d statistics and the l ower energy thresho ld that 



is now available at the Pierre Auger Observatory (jMathes et al. 20111 : ISanchez et al. 20111 ) . This 
will provide further constraints helping to understand the origin of cosmic rays in the energy range 
0.1 < E/EeV < 10. 
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Appendix A : Large scale anisotropies in local coordinates 
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Fig. 12. — Effect of large scale anisotropies in local coordinates (left : as a function of sin 2 0, right : as a 
function of ip) for an observer located at the Earth latitude 4itc = —35.2° of the Pierre Auger Observatory. 



To study the angular distribution in local coordinates for different anisotropic angular distribu- 
tions &(a, 5) in celestial coordinates, we restrict ourselves, without loss of generalities, to the case 
of full detection efficiency (e(0,<p,E) = 1). Then, the instantaneous arrival direction distribution 
in local coordinates reads : 

d 3 7V 

ocsinfl cosfl g(fl,p,a°). (29) 

dudcpdar 

&(6,ip,a°) is the underlying angular distribution of cosmic rays, expressed in local coordinates. 
In case of isotropy, <1> is constant so that once integrated over ip and cv°, the arrival direction 
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distribution is such that diV/dsin 2 9 is also constant. On the other hand, in case of a dipolar 
distribution for instance, $ is proportional to l + rd(0, tp, a ) • n(9, tp), where n is here a unit vector 
in local coordinates, and d the dipole unit vector pointing towards (ad, 5d) and expressed in local 
coordinates by means of Eqn.|H To quantify the distortions induced by a dipole in the d-ZV/dsin 2 9 
distribution, we define A(diV/dsin 2 9) such that : 

A I a at I A -2m 1 ( djVgggZe/d sin 2 9 - dN lso /d sin 2 9 \ 

A(diV/dsm 9) = - ^ diVWdWg J" (30) 

Once multiplied by the dipole amplitude r, A(dA r /dsin 2 0) gives directly the relative changes in the 
diV/dsin 2 # distribution with respect to isotropy. Carrying out integrations over tp and a yields 
to : 

A(diV/dsin 2 60 = N °fP° le sin4it e sin<5 d cosfl, (31) 

Nq^so 

where both intensity normalisations Nq^o and NQ^ipole are tuned to guarantee the same number 
of events observed in the covered region of the sky for each underlying angular distribution. This 
result is shown in the left panel of Fig. [T2l for the latitude £ s it e = —35.2° of the Pierre Auger 
Observatory and for different dipole directions. Within the zenithal range [0°,55°] considered in 
this article, the relative changes - maximal for 5^ = ±90° - amount at most to ~ ±15%. So, even 
for an amplitude r as large as 10%, the relative changes in diV/dsin 2 9 would be within ~ ±1.5%, 
variation which - given the available statistics - is sufficiently low to be considered as negligible. 
Besides, the same calculation applied to the case of a symmetric quadrupolar anisotropy shows 
that the variation of A(diV/dsin 2 9) is less than ~ 0.1%, thus being negligible. Consequently, the 
distribution in d-/V/dsin 2 9 can be considered at first order as insensitive to large scale anisotropics, 
so that any significant deviation from a uniform distribution provides an empirical measurement of 
the zenithal dependence of the detection efficiency 

It is worth noting that the azimuthal distribution averaged over time is, on the other hand, 
sensitive to large scale anisotropics. Repeating the same calculation and integrating now over 9 (in 
this example between and 60°) and a yields the A(dN/dtp) relative changes : 

A I A AT I A \ N 0,dipole SU1 5 d COS 4ite „ . „ fX . \ , QrA 

A(dA7dv3) = ~jrz — 7tan4ite + WSsm tp . (32) 

This function is shown in the right panel of Fig. [T2l for 5^ = 90° (dashed line) and 8^ = —90° (dotted 
line). The amplitude of the dipole wave is now ~ 0.5. As well, the influence of a quadrupole on 
A(dA r /d(/?) is illustrated by the dashed-dotted line (oblate symmetric quadrupole in this example). 
Since, at the Earth latitude of the Pierre Auger Observatory, any genuine large scale pattern which 
depends on the declination translates into azimuthal modulations of the event rate similar to the 
ones induced by experimental effects, it is thus mandatory to model accurately the dependence on 
azimuth of the detection efficiency for disentangling local from celestial effects. 



-37- 



Appendix B : Modulation of the detection efficiency induced by a tilted array 

To estimate the modulation of the detection efficiency induced by a tilted array, we consider 
here that in the absence of tilt, the corresponding detection efficiency function e no tiit depends only 
on the energy and the zenith angle and can be parameterised in a good approximation as : 

E 3 

e w m(E,e) = E , + EUe y (33) 
Eqs{0) is the zenithal-dependent energy at which c no tilt 

(E, 6) = 0.5. In case of a tilted array, this 
parameter depends also on the azimuth angle, which is then the source of the azimuthal modulation 
of the detection efficiency. To understand this, it is useful to consider for any given shower with 
parameters (E, 0, ip) the circle in the shower plane corresponding to the region in which a signal S 
larger than some specified threshold value So is expected. Let ro(C) denote the radius of this circle, 
£ being the tilt angle of the SD array. The detection efficiency, and hence also the parameter Eq^, 
is ultimately a function of the average number of detectors contained in the projection of this circle 
into the ground, given by : 

2 

(n det ) (5 > So) oc r , (34) 

n, z |n_|_ • n| 

where h = 1.5 km is the nominal separation between surface detectors. The radii ro(C) obtained 
with the tilted array leading to the same value of (ndet) can be related to ro(C = 0) through : 

tfO-*<-0)^. (35) 
Hence, we can obtain the relation between the energies 

£o.5 with tilt (E$%) and without tilt (E , 5 ) 
by comparing the cosmic ray energies required to get the value So a t radius ro(C) and at radius 
r o(C = 0)- Approximating the lateral distribution function of the signal near the radius r$ as a 
power law S(r) oc Er~ 3 , we obtain the following relation : 

E$l(0,<p)=E o . s (6)( ?}® n . ) ^£ O .5(#)[1 + Ctan0cos(^-^ O )] 3 . (36) 



MC = o) 

Then, subtracting e no tnt to e t ut leads to Eqn. [9l 



Appendix C : Determination of upper limits on dipole amplitudes 



To determine upper limits on the dipole amplitud es, Linsley de scribed the procedure to follow 
in the case of first harmonic analysis in right ascension (jLinslev 19751 ) . We adapt here this procedure 
to the case of the dipolar reconstruction adopted in section 15.21 



Here, the data set is supposed to have been drawn at random from an underlying dipolar 
distribution characterised by d, whose value is unknown. In the limit of large number of events, 
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the joint p.d.f. PD x ,D Y ,D z (d x ,d y ,d z ) can be factorised in terms of three Gaussian distributions 
N(di - di,ai) : 



PD x ,D Y ,D z {d x ,dy,d z ;d x ,dy,d z ) = N(d x - d x ,a)N(d y - d y ,a)N(d z - d z ,a z ). 



(37) 



The joint p.d.f. PR.A,A(r, S, a) expressing the dipole components in spherical coordinates is then 
obtained by performing the Jacobian transformation : 



d(d x ,d y ,d ZJ 



d(r, 8, a) 
r 2 cos 8 
(2ir) 3 / 2 a 2 a : 



PD x ,D Y ,D z (d x (r, 5, a),d y (r, 8, a), d z (r, 8, a)) 
(r sin 8 — d sin 5d) 2 



cxp 



2a? 



x exp 
x exp 



(r cos 8 cos a — d cos Sd cos ay) 

(r cos 8 sin a — d cos 6 a sin a^Y 
2a 2 



(38) 



Each analysed data set having been selected at random from an ensemble in which all possible 
values of d are equally represented, the various d, 8^ and combinations have relative probability 
PR,A,A(r, 8,a; d, 8d, ad) /pr,a,a(j, 8 ,a; d = 0). This allows us to define the joint p.d.f. pr,a,a by 
requiring this ratio to be normalised to unity : 



PR,A,A(r, 5, a; d, 8 d ,a d ) 



K(r, 8) exp 



rd cos 8 cos 8d cos (a — ay) 



(7- 



x exp 



rd sin 8 sin 8d d 2 cos 2 8d d 2 sin 2 



2a 2 



2al 



(39) 



where the normalisation reads : 



K(r, 8) 



2ir In 



rd cos 8 cos 8 n 



o- 



dd d8d exp 



d 2, cos^ 8d d 2 sin 2 8d rd sin 8 sin 8d 



2a 2 



2o1 



+ 



(40) 



In is here the modified Bessel function of the first kind with order 0. Integration of Pr,a,a over 8d 
and ad yields the pr p.d.f., from which upper limits on d can be obtained within a confidence level 
C.L. by inverting the relation : 



dfp R (r,8;d UL ) = C.L. 



(41) 



T data 



Due to the non-uniform directional exposure in declination, the resulting upper limits actually 
depend on the declination through the dependence of Pr on 8. In practice, this dependence is 
small, which is why we presented in section [7] upper limits averaged over the declination. 
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Appendix D : Determination of upper limits on quadrupole amplitudes 

To determine upper limits on quadrupole amplitudes, we rely on Monte-Carlo simulations. For 
each possible amplitude A + (/3), we estimate the p.d.f. j>a + (A + ; A+) (pb(/3; P)) with a given number 
of events N and a given exposure 60. The amplitude \% L such that fj° dA+ p A (A+; A^ L ) = C.L. 
is a relevant upper limit (and respectively for /3 ). 

Alternatively to the previous procedure used to derive upper limits on dipole amplitudes, this 
procedure can lead to upper limits tighter than the upper bounds for isotropy A +i gg when the 
measured values of \+ t data are smaller than the expected average for isotropy. To cope with this 
undesired behaviour, the upper limits presented in section [7] are defined as max(A +i 99, A+ L ). 

REFERENCES 

R. U. Abbasi et al. (The IceCube Collaboration), ApJ 718 (2010) 194 

R. U. Abbasi et al. (The HiRes Collaboration), Phys. Rev. Lett. 104 (2010) 161101 

A. A. Abdo et al. (The Milagro Collaboration), ApJ 698 (2009) 2121 

M. Aglietta et al. (The EAS-TOP Collaboration), ApJL 692 (2009) L130-L133 

M. Amenomori et al. (The Tibet AS7 Collaboration), Science 314 (2006) 439 

R. Beck, Space Sci. Rev. 99 (2001) 243 

V. S. Berezinsky, S. I. Grigorieva, B. I. Hnatyk, Astropart. Phys. 21 (2004) 617625 
V. S. Berezinsky, A. Z. Gazizov, S. I. Grigorieva, Phys. Rev. D 74 (2006) 043005 
P. Billoir, O. Deligny, JCAP 02 (2008) 009 

D. J. Bird et al, (The Fly's Eye Collaboration), Phys. Rev. Lett. 71 (1993) 3401 
G. R. Blumenthal, Phys. Rev. Dl (1970) 1596 

C. Bonifazi, A. Letessier-Selvon, E.M. Santos, Astropart. Phys. 28 (2008) 523 

C. Bonifazi for the Pierre Auger Collaboration, Nucl. Phys. Proc. Suppl. 190 (2009) 20 
A. Calvez, A. Kusenko, S. Nagataki, Phys. Rev. Lett. 105 (2010) 091101 

J. Candia, S. Mollerach, E. Roulet, JCAP 0305 (2003) 003 
A. H. Compton, I. A. Getting, Phys. Rev. 47 (1935) 817 

D. Cutler, D. Groom, Nature 322 (1986) 434 



-40- 



D. Eichler, M. Pohl, ApJ 742 (2011) 114 

R. Engel for the Pierre Auger Collaboration, Proceedings of the 30th ICRC, Merida, 2007 
G. Farrar, R. Jansson, ApJ 757 (2012) 14 

G. Giacinti et al, JCAP 07 (2012) 031 

D. Harari, S. Mollerach, E. Roulet, JCAP 11 (2010) 033 
J. Hersil et al, Phys. Rev. Lett. 6 (1961) 22 
A. M. Hillas, Phys. Lett. 24A (1967) 677 

C. C. H. Jui et al (The Telescope Array Collaboration), Proceedings of the APS meeting, 2011, 
larXrv:1110.0133l 

M. Kachelriess, P. Serpico, Phys. Lett. B 640 (2006) 225-229 

M. A. Lawrence, R. J. O. Reid, A. A. Watson, J. Phys. G 17 (1991) 733 

T.-P. Li, Y.-Q. Ma, ApJ 272 (1983) 317 

J. Linsley, Proceedings of the 8th ICRC, Jaipur, vol. 4, 1963, p. 77 
J. Linsley, Phys. Rev. Lett. 34 (1975) 221101 

H. J. Mathes for the Pierre Auger Collaboration, Proceedings of the 32nd ICRC, Beijing, 2011 
M. Nagano et al, J. Phys. G 18 (1992) 423 

R. Pesce for the Pierre Auger Collaboration, Proceedings of the 32nd ICRC, Beijing, 2011 

The Pierre Auger Collaboration, Nucl. Instr. and Meth. A 523 (2004) 50 

The Pierre Auger Collaboration, Phys. Rev. Lett. 101 (2008) 061101 

The Pierre Auger Collaboration, Astropart. Phys. 32 (2009) 89 

The Pierre Auger Collaboration, Phys. Lett. B 685 (2010a) 239 

The Pierre Auger Collaboration, Nucl. Instr. and Meth. A 613 (2010b) 29 

The Pierre Auger Collaboration, Phys. Rev. Lett. 104 (2010c) 091101 

The Pierre Auger Collaboration, Astropart. Phys. 34 (2011a) 627-639 

The Pierre Auger Collaboration, JCAP 11 (2011b) 022 

The Pierre Auger Collaboration, in preparation (2012a) 



-41 - 



The Pierre Auger Collaboration, submitted to ApJ (2012b) 

M. S. Pshirkov et al, ApJ 738 (2011) 192 

V. Ptuskin et al, Astron. Astrophys. 268 (1993) 726 

F. Sanchez for the Pierre Auger Collaboration, Proceedings of the 32nd ICRC, Beijing, 2011 

P. Sommers, Astropart. Phys. 14 (2001) 71 

K. O. Thielheim, W. Langhoff, J. Phys. A (1968) 694 

V. N. Zirakashvili et al, AstL. 24 (1998) 139 



This preprint was prepared with the AAS IATjrjX macros v5.2. 



